{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Transforms and Resampling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook explains how to apply transforms to images, and how to perform image resampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "import SimpleITK as sitk\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "from ipywidgets import interact, fixed\n",
    "\n",
    "# Utility method that either downloads data from the Girder repository or\n",
    "# if already downloaded returns the file name for reading from disk (cached data).\n",
    "%run update_path_to_download_script\n",
    "from downloaddata import fetch_data as fdata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating and Manipulating Transforms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A number of different spatial transforms are available in SimpleITK.\n",
    "\n",
    "The simplest is the Identity Transform.  This transform simply returns input points unaltered."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dimension = 2\n",
    "\n",
    "print('*Identity Transform*')\n",
    "identity = sitk.Transform(dimension, sitk.sitkIdentity)\n",
    "print('Dimension: ' + str(identity.GetDimension()))\n",
    "\n",
    "# Points are always defined in physical space\n",
    "point = (1.0, 1.0)\n",
    "def transform_point(transform, point):\n",
    "    transformed_point = transform.TransformPoint(point)\n",
    "    print('Point ' + str(point) + ' transformed is ' + str(transformed_point))\n",
    "\n",
    "transform_point(identity, point)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Transform are defined by two sets of parameters, the *Parameters* and *FixedParameters*.  *FixedParameters* are not changed during the optimization process when performing registration.  For the TranslationTransform, the Parameters are the values of the translation Offset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('*Translation Transform*')\n",
    "translation = sitk.TranslationTransform(dimension)\n",
    "\n",
    "print('Parameters: ' + str(translation.GetParameters()))\n",
    "print('Offset:     ' + str(translation.GetOffset()))\n",
    "print('FixedParameters: ' + str(translation.GetFixedParameters()))\n",
    "transform_point(translation, point)\n",
    "\n",
    "print('')\n",
    "translation.SetParameters((3.1, 4.4))\n",
    "print('Parameters: ' + str(translation.GetParameters()))\n",
    "transform_point(translation, point)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The affine transform is capable of representing translations, rotations, shearing, and scaling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('*Affine Transform*')\n",
    "affine = sitk.AffineTransform(dimension)\n",
    "\n",
    "print('Parameters: ' + str(affine.GetParameters()))\n",
    "print('FixedParameters: ' + str(affine.GetFixedParameters()))\n",
    "transform_point(affine, point)\n",
    "\n",
    "print('')\n",
    "affine.SetTranslation((3.1, 4.4))\n",
    "print('Parameters: ' + str(affine.GetParameters()))\n",
    "transform_point(affine, point)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A number of other transforms exist to represent non-affine deformations, well-behaved rotation in 3D, etc. See the [Transforms](22_Transforms.ipynb) tutorial for more information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Applying Transforms to Images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a function to display the images that is aware of image spacing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def myshow(img, title=None, margin=0.05, dpi=80):\n",
    "    nda = sitk.GetArrayViewFromImage(img)\n",
    "    spacing = img.GetSpacing()\n",
    "        \n",
    "    ysize = nda.shape[0]\n",
    "    xsize = nda.shape[1]\n",
    "      \n",
    "    figsize = (1 + margin) * ysize / dpi, (1 + margin) * xsize / dpi\n",
    "\n",
    "    fig = plt.figure(title, figsize=figsize, dpi=dpi)\n",
    "    ax = fig.add_axes([margin, margin, 1 - 2*margin, 1 - 2*margin])\n",
    "    \n",
    "    extent = (0, xsize*spacing[1], 0, ysize*spacing[0])\n",
    "    \n",
    "    t = ax.imshow(nda,\n",
    "            extent=extent,\n",
    "            interpolation='hamming',\n",
    "            cmap='gray',\n",
    "            origin='lower')\n",
    "    \n",
    "    if(title):\n",
    "        plt.title(title)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a grid image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = sitk.GridSource(outputPixelType=sitk.sitkUInt16,\n",
    "    size=(250, 250),\n",
    "    sigma=(0.5, 0.5),\n",
    "    gridSpacing=(5.0, 5.0),\n",
    "    gridOffset=(0.0, 0.0),\n",
    "    spacing=(0.2,0.2))\n",
    "myshow(grid, 'Grid Input')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To apply the transform, a resampling operation is required."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def resample(image, transform):\n",
    "    # Output image Origin, Spacing, Size, Direction are taken from the reference\n",
    "    # image in this call to Resample\n",
    "    reference_image = image\n",
    "    interpolator = sitk.sitkCosineWindowedSinc\n",
    "    default_value = 100.0\n",
    "    return sitk.Resample(image, reference_image, transform,\n",
    "                         interpolator, default_value)\n",
    "\n",
    "translation.SetOffset((3.1, 4.6))\n",
    "transform_point(translation, point)\n",
    "resampled = resample(grid, translation)\n",
    "myshow(resampled, 'Resampled Translation')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happened?  The translation is positive in both directions.  Why does the output image move down and to the left?  It important to keep in mind that a transform in a resampling operation defines *the transform from the output space to the input space*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "translation.SetOffset(-1*np.array(translation.GetParameters()))\n",
    "transform_point(translation, point)\n",
    "resampled = resample(grid, translation)\n",
    "myshow(resampled, 'Inverse Resampled')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An affine (line preserving) transformation, can perform translation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def affine_translate(transform, x_translation=3.1, y_translation=4.6):\n",
    "    new_transform = sitk.AffineTransform(transform)\n",
    "    new_transform.SetTranslation((x_translation, y_translation))\n",
    "    resampled = resample(grid, new_transform)\n",
    "    myshow(resampled, 'Translated')\n",
    "    return new_transform\n",
    "    \n",
    "affine = sitk.AffineTransform(dimension)\n",
    "\n",
    "interact(affine_translate, transform=fixed(affine), x_translation=(-5.0, 5.0), y_translation=(-5.0, 5.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or scaling:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def affine_scale(transform, x_scale=3.0, y_scale=0.7):\n",
    "    new_transform = sitk.AffineTransform(transform)\n",
    "    matrix = np.array(transform.GetMatrix()).reshape((dimension,dimension))\n",
    "    matrix[0,0] = x_scale\n",
    "    matrix[1,1] = y_scale\n",
    "    new_transform.SetMatrix(matrix.ravel())\n",
    "    resampled = resample(grid, new_transform)\n",
    "    myshow(resampled, 'Scaled')\n",
    "    print(matrix)\n",
    "    return new_transform\n",
    "\n",
    "affine = sitk.AffineTransform(dimension)\n",
    "\n",
    "interact(affine_scale, transform=fixed(affine), x_scale=(0.2, 5.0), y_scale=(0.2, 5.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or rotation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def affine_rotate(transform, degrees=15.0):\n",
    "    parameters = np.array(transform.GetParameters())\n",
    "    new_transform = sitk.AffineTransform(transform)\n",
    "    matrix = np.array(transform.GetMatrix()).reshape((dimension,dimension))\n",
    "    radians = -np.pi * degrees / 180.\n",
    "    rotation = np.array([[np.cos(radians), -np.sin(radians)],[np.sin(radians), np.cos(radians)]])\n",
    "    new_matrix = np.dot(rotation, matrix)\n",
    "    new_transform.SetMatrix(new_matrix.ravel())\n",
    "    resampled = resample(grid, new_transform)\n",
    "    print(new_matrix)\n",
    "    myshow(resampled, 'Rotated')\n",
    "    return new_transform\n",
    "    \n",
    "affine = sitk.AffineTransform(dimension)\n",
    "\n",
    "interact(affine_rotate, transform=fixed(affine), degrees=(-90.0, 90.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or shearing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def affine_shear(transform, x_shear=0.3, y_shear=0.1):\n",
    "    new_transform = sitk.AffineTransform(transform)\n",
    "    matrix = np.array(transform.GetMatrix()).reshape((dimension,dimension))\n",
    "    matrix[0,1] = -x_shear\n",
    "    matrix[1,0] = -y_shear\n",
    "    new_transform.SetMatrix(matrix.ravel())\n",
    "    resampled = resample(grid, new_transform)\n",
    "    myshow(resampled, 'Sheared')\n",
    "    print(matrix)\n",
    "    return new_transform\n",
    "\n",
    "affine = sitk.AffineTransform(dimension)\n",
    "\n",
    "interact(affine_shear, transform=fixed(affine), x_shear=(0.1, 2.0), y_shear=(0.1, 2.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Composite Transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is possible to compose multiple transform together into a single transform object.  With a composite transform, multiple resampling operations are prevented, so interpolation errors are not accumulated.  For example, an affine transformation that consists of a translation and rotation,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "translate = (8.0, 16.0)\n",
    "rotate = 20.0\n",
    "\n",
    "affine = sitk.AffineTransform(dimension)\n",
    "affine = affine_translate(affine, translate[0], translate[1])\n",
    "affine = affine_rotate(affine, rotate)\n",
    "\n",
    "resampled = resample(grid, affine)\n",
    "myshow(resampled, 'Single Transform')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "can also be represented with two Transform objects applied in sequence with a Composite Transform,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "composite = sitk.Transform(dimension, sitk.sitkComposite)\n",
    "translation = sitk.TranslationTransform(dimension)\n",
    "translation.SetOffset(-1*np.array(translate))\n",
    "composite.AddTransform(translation)\n",
    "affine = sitk.AffineTransform(dimension)\n",
    "affine = affine_rotate(affine, rotate)\n",
    "\n",
    "composite.AddTransform(translation)\n",
    "composite = sitk.Transform(dimension, sitk.sitkComposite)\n",
    "composite.AddTransform(affine)\n",
    "\n",
    "resampled = resample(grid, composite)\n",
    "myshow(resampled, 'Two Transforms')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Beware*, tranforms are non-commutative -- order matters!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "composite = sitk.Transform(dimension, sitk.sitkComposite)\n",
    "composite.AddTransform(affine)\n",
    "composite.AddTransform(translation)\n",
    "\n",
    "resampled = resample(grid, composite)\n",
    "myshow(resampled, 'Composite transform in reverse order')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Resampling\n",
    "\n",
    "<img src=\"resampling.svg\"/><br><br>\n",
    "\n",
    "Resampling as the verb implies is the action of sampling an image, which itself is a sampling of an original continuous signal.\n",
    "\n",
    "Generally speaking, resampling in SimpleITK involves four components:\n",
    "1. Image - the image we resample, given in coordinate system $m$.\n",
    "2. Resampling grid - a regular grid of points given in coordinate system $f$ which will be mapped to coordinate system $m$.\n",
    "2. Transformation $T_f^m$ - maps points from coordinate system $f$ to coordinate system $m$, $^mp = T_f^m(^fp)$.\n",
    "3. Interpolator - method for obtaining the intensity values at arbitrary points in coordinate system $m$ from the values of the points defined by the Image.\n",
    "\n",
    "\n",
    "While SimpleITK provides a large number of interpolation methods, the two most commonly used are ```sitkLinear``` and ```sitkNearestNeighbor```. The former is used for most interpolation tasks, a compromise between accuracy and computational efficiency. The later is used to interpolate labeled images representing a segmentation, it is the only interpolation approach which will not introduce new labels into the result.\n",
    "\n",
    "SimpleITK's procedural API provides three methods for performing resampling, with the difference being the way you specify the resampling grid:\n",
    "\n",
    "1. ```Resample(const Image &image1, Transform transform, InterpolatorEnum interpolator, double defaultPixelValue, PixelIDValueEnum outputPixelType)```\n",
    "2. ```Resample(const Image &image1, const Image &referenceImage, Transform transform, InterpolatorEnum interpolator, double defaultPixelValue, PixelIDValueEnum outputPixelType)```\n",
    "3. ```Resample(const Image &image1, std::vector< uint32_t > size, Transform transform, InterpolatorEnum interpolator, std::vector< double > outputOrigin, std::vector< double > outputSpacing, std::vector< double > outputDirection, double defaultPixelValue, PixelIDValueEnum outputPixelType)```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def resample_display(image, euler2d_transform, tx, ty, theta):\n",
    "    euler2d_transform.SetTranslation((tx, ty))\n",
    "    euler2d_transform.SetAngle(theta)\n",
    "    \n",
    "    resampled_image = sitk.Resample(image, euler2d_transform)\n",
    "    plt.imshow(sitk.GetArrayFromImage(resampled_image))\n",
    "    plt.axis('off')    \n",
    "    plt.show()\n",
    "\n",
    "logo = sitk.ReadImage(fdata('SimpleITK.jpg'))\n",
    "\n",
    "euler2d = sitk.Euler2DTransform()\n",
    "# Why do we set the center?\n",
    "euler2d.SetCenter(logo.TransformContinuousIndexToPhysicalPoint(np.array(logo.GetSize())/2.0))\n",
    "interact(resample_display, image=fixed(logo), euler2d_transform=fixed(euler2d), tx=(-128.0, 128.0,2.5), ty=(-64.0, 64.0), theta=(-np.pi/4.0,np.pi/4.0));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Common Errors\n",
    "\n",
    "It is not uncommon to end up with an empty (all black) image after resampling. This is due to:\n",
    "1. Using wrong settings for the resampling grid, not too common, but does happen.\n",
    "2. Using the inverse of the transformation $T_f^m$. This is a relatively common error, which is readily addressed by invoking the transformations ```GetInverse``` method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining the Resampling Grid\n",
    "\n",
    "In the example above we arbitrarily used the original image grid as the resampling grid. As a result, for many of the transformations the resulting image contained black pixels, pixels which were mapped outside the spatial domain of the original image and a partial view of the original image.\n",
    "\n",
    "If we want the resulting image to contain all of the original image no matter the transformation, we will need to define the resampling grid using our knowledge of the original image's spatial domain and the **inverse** of the given transformation. \n",
    "\n",
    "Computing the bounds of the resampling grid when dealing with an affine transformation is straightforward. An affine transformation preserves convexity with extreme points mapped to extreme points. Thus we only need to apply the **inverse** transformation to the corners of the original image to obtain the bounds of the resampling grid.\n",
    "\n",
    "Computing the bounds of the resampling grid when dealing with a BSplineTransform or DisplacementFieldTransform is more involved as we are not guaranteed that extreme points are mapped to extreme points. This requires that we apply the **inverse** transformation to all points in the original image to obtain the bounds of the resampling grid.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "euler2d = sitk.Euler2DTransform()\n",
    "# Why do we set the center?\n",
    "euler2d.SetCenter(logo.TransformContinuousIndexToPhysicalPoint(np.array(logo.GetSize())/2.0))\n",
    "\n",
    "tx = 64\n",
    "ty = 32\n",
    "euler2d.SetTranslation((tx, ty))\n",
    "\n",
    "extreme_points = [logo.TransformIndexToPhysicalPoint((0,0)), \n",
    "                  logo.TransformIndexToPhysicalPoint((logo.GetWidth(),0)),\n",
    "                  logo.TransformIndexToPhysicalPoint((logo.GetWidth(),logo.GetHeight())),\n",
    "                  logo.TransformIndexToPhysicalPoint((0,logo.GetHeight()))]\n",
    "inv_euler2d = euler2d.GetInverse()\n",
    "\n",
    "extreme_points_transformed = [inv_euler2d.TransformPoint(pnt) for pnt in extreme_points]\n",
    "min_x = min(extreme_points_transformed)[0]\n",
    "min_y = min(extreme_points_transformed, key=lambda p: p[1])[1]\n",
    "max_x = max(extreme_points_transformed)[0]\n",
    "max_y = max(extreme_points_transformed, key=lambda p: p[1])[1]\n",
    "\n",
    "# Use the original spacing (arbitrary decision).\n",
    "output_spacing = logo.GetSpacing()\n",
    "# Identity cosine matrix (arbitrary decision).   \n",
    "output_direction = [1.0, 0.0, 0.0, 1.0]\n",
    "# Minimal x,y coordinates are the new origin.\n",
    "output_origin = [min_x, min_y]\n",
    "# Compute grid size based on the physical size and spacing.\n",
    "output_size = [int((max_x-min_x)/output_spacing[0]), int((max_y-min_y)/output_spacing[1])]\n",
    "\n",
    "resampled_image = sitk.Resample(logo, output_size, euler2d, sitk.sitkLinear, output_origin, output_spacing, output_direction)\n",
    "plt.imshow(sitk.GetArrayViewFromImage(resampled_image))\n",
    "plt.axis('off')    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Are you puzzled by the result? Is the output just a copy of the input? Add a rotation to the code above and see what happens (```euler2d.SetAngle(0.79)```)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Resampling at a set of locations \n",
    "\n",
    "In some cases you may be interested in obtaining the intensity values at a set of points (e.g. coloring the vertices of a mesh model segmented from an image).\n",
    "\n",
    "The code below generates a random point set in the image and resamples the intensity values at these locations. It is written so that it works for all image-dimensions and types (scalar or vector pixels)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = logo\n",
    "\n",
    "# Generate random samples inside the image, we will obtain the intensity/color values at these points.\n",
    "num_samples = 10\n",
    "physical_points = []\n",
    "for pnt in zip(*[list(np.random.random(num_samples)*sz) for sz in img.GetSize()]):\n",
    "    physical_points.append(img.TransformContinuousIndexToPhysicalPoint(pnt))\n",
    "\n",
    "# Create an image of size [num_samples,1...1], actual size is dependent on the image dimensionality. The pixel\n",
    "# type is irrelevant, as the image is just defining the interpolation grid (sitkUInt8 has minimal memory footprint).\n",
    "interp_grid_img = sitk.Image([num_samples] + [1]*(img.GetDimension()-1), sitk.sitkUInt8)\n",
    "\n",
    "# Define the displacement field transformation, maps the points in the interp_grid_img to the points in the actual\n",
    "# image.\n",
    "displacement_img = sitk.Image([num_samples] + [1]*(img.GetDimension()-1), sitk.sitkVectorFloat64, img.GetDimension())\n",
    "for i, pnt in enumerate(physical_points):\n",
    "     displacement_img[[i] + [0]*(img.GetDimension()-1)] = np.array(pnt) - np.array(interp_grid_img.TransformIndexToPhysicalPoint([i] + [0]*(img.GetDimension()-1)))\n",
    "\n",
    "# Actually perform the resampling. The only relevant choice here is the interpolator. The default_output_pixel_value\n",
    "# is set to 0.0, but the resampling should never use it because we expect all points to be inside the image and this\n",
    "# value is only used if the point is outside the image extent.\n",
    "interpolator_enum = sitk.sitkLinear\n",
    "default_output_pixel_value = 0.0\n",
    "output_pixel_type = sitk.sitkFloat32 if img.GetNumberOfComponentsPerPixel()==1 else sitk.sitkVectorFloat32\n",
    "resampled_points = sitk.Resample(img, interp_grid_img, sitk.DisplacementFieldTransform(displacement_img), \n",
    "                                 interpolator_enum, default_output_pixel_value, output_pixel_type)\n",
    "\n",
    "# Print the interpolated values per point\n",
    "for i in range(resampled_points.GetWidth()):\n",
    "      print(str(physical_points[i]) + ': ' + str(resampled_points[[i] + [0]*(img.GetDimension()-1)]) + '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=\"red\">Homework:</font> creating a color mesh\n",
    "\n",
    "You will now use the code for resampling at arbitrary locations to create a colored mesh.\n",
    "\n",
    "Using the color image of the [visible human](https://en.wikipedia.org/wiki/Visible_Human_Project) head [`img = sitk.ReadImage(fdata('vm_head_rgb.mha'))`]:\n",
    "1. Implement the [marching cubes algorithm](https://en.wikipedia.org/wiki/Marching_cubes) to obtain the set of triangles corresponding to the iso-surface of structures of interest (skin, white matter,...). The original paper is available from the authors [here](http://marchingcubes.org/images/f/f9/MarchingCubes.pdf).\n",
    "2. Find the color associated with each of the triangle vertices using the code above.\n",
    "3. Save the data using the ASCII version of the [PLY](https://en.wikipedia.org/wiki/PLY_(file_format)), Polygon File Format (a.k.a. Stanford Triangle Format).\n",
    "4. Use [meshlab](http://www.meshlab.net/) to view your creation."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
